上一篇我們提到了二元分類的logistic regression，今天來看一下多類別的實做方式。

補充一下昨天沒有寫清楚的部份

這邊的Yn是，第n筆資料是第一個類別的機率，因為是二元分類，所以(1 - Yn)就是第二個類別的機率。

那多元分類有幾個跟二元稍微不同的地方，首先二元分類的Yn是以sigmoid function，而在多元分類中，會以softmax來表示，寫成

意思就是第n筆資料是第k個類別的機率是多少，其中

知道這個之後就可以簡單的寫下likelihood

並且得到error function

而昨天提到的牛頓法，我們需要的就是這個error function的一次與二次微分

接著我們就可以利用這兩個式子來寫程式囉！我用的資料集在這。

import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt
from math import isnan

data =  np.loadtxt(open('train.csv' , 'r'), delimiter = ',')
target = data[:,0:3]
feature = data[:,3:]

def phi(x):
        x = np.reshape(x,(len(x),1))
        return x
        
def y(n,k,w,X): #softmax
        s = np.float64(0.)
	ak = w[k].T.dot(phi(X[n]))
        for j in range(3):
                aj = w[j].T.dot(phi(X[n]))
                s += np.nan_to_num(np.exp(aj - ak))
        s = np.nan_to_num(s)
        return 1./s
        
def gradient(w,k,t,X):
	output = np.zeros((len(w[0]),1))
	for n in range(len(X)):
		scale = y(n,k,w,X) - t[:,k][n] #Ynk - Tnk
		output += scale * phi(X[n])
	return output
    
def hessian(w,k,X):
	output = np.zeros((len(w[0]),len(w[0])))
	for n in range(len(X)): 
		scale = y(n,k,w,X) * (1 - y(n,k,w,X))
		output += scale * (phi(X[n]).dot(phi(X[n]).T)) 
	return output
    
def error(w,t,X):
	s = np.float64(0.)
	for n in range(len(X)):
		for k in range(3):
			if t[:,k][n] != 0.:
				s += np.nan_to_num(np.log(y(n,k,w,X)))
	return -1*s
    
def classify(w,x):
	softmax = []
	for k in range(3):
		s = np.float64(0.)
		ak = w[k].T.dot(phi(x))
		for j in range(3):
			aj = w[j].T.dot(phi(x))
			s += np.nan_to_num(np.exp(aj - ak))
		softmax += [1./s]
	return softmax.index(max(softmax))

w = np.zeros((3,len(phi(feature[0])),1))
cee = []
acc = []
while True:
	e = error(w,target, feature)
	acc += [accuracy(w,target,feature)]
	cee += [np.reshape(e,1)]
	if e < 0.001:
		break
	for k in range(3):
		w[k] = w[k] - inv(hessian(w,k,feature)).dot(gradient(w,k,target,feature))

epoch = [i for i in range(len(cee))]
plt.xlabel('epoch')
plt.ylabel('cross entropy error')
plt.plot(epoch, cee, linestyle = '-')
plt.xticks(epoch)
plt.savefig('cross_entropy_error.png')

結果大概會得到這樣的圖，確定cross entropy error的確有在下降，也證明這個演算法是有用的！